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Abstract 

A theoretical study has been conducted on the effects of cloud horizontal inhomogeneity on 
cloud albedo bias. A two-dimensional (2D) version of the Spherical Harmonic Discrete Ordinate 
Method (SHDOM) is used to estimate the albedo bias of the plane parallel (PP-IPA) and 
independent pixel (IPA-2D) approximations for a wide range of 2D cloud fields obtained from 
LANDSAT. They include single layer trade cumulus, open and closed cell broken stratocumulus, 
and solid stratocumulus boundary layer cloud fields over ocean. Findings are presented on a 
variety of averaging scales and are summarized as a function of cloud fraction, mean cloud 
optical depth, cloud aspect ratio, standard deviation of optical depth, and the gamma function 
parameter v (a measure of the width of the optical depth distribution). Biases are found to be 
small for small cloud fraction or mean optical depth, where the cloud fields under study behave 
linearly. They are large (up to 0.20 for PP-IPA bias, -0.12 for IPA-2D bias) for large v. On a 
scene average basis PP-IPA bias can reach 0.30, while IPA-2D bias reaches its largest magnitude 
at -0.07. Biases due to horizontal transport (IPA-2D) are much smaller than PP-IPA biases but 
account for 20% RMS of the bias overall. 

Limitations of this work include the particular cloud field set used, assumptions of conservative 
scattering, constant cloud droplet size, no gas absorption or surface reflectance, and restriction to 
2D radiative transport. The LANDSAT data used may also be affected by radiative smoothing. 

1. Introduction 

Recent results by Cahalan et al. (1994a, b) have shown that the albedo bias of the independent pixel approx- 
imation (IPA-3D) for marine boundary layer clouds is about 1% of the typical albedo, while the plane parallel 
(PP-IPA) bias is on the order of 15%. That study included only overcast conditions at one location off the coast of 
California. This note extends Cahalan’s work to a wide range of marine boundary layer cloud conditions, including 
broken stratus and trade cumulus cases, found in 45 LANDSAT scenes of clouds over ocean, ranging from about 
20 S to 40 N latitude and from 15 to 150 W longitude. Computations of the radiance field of these scenes using 
both IPA and 2D assumptions are performed at 0.83 //m. These calculations allow estimation of the error in the 
PP and IP approximations for a wide range of cloud parameters, complementing Cahalan’s results for overcast 
marine stratocumulus. 
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A limitation of the current work is that the IPA is used initially to derive the LANDSAT optical depth distri- 
butions which form the reference field for later testing of the IPA and 2D calculations. Furthermore, since only 2D 
cloud fields are considered here (effectively, cloud streets), the full effect of realistic cloud fields has not been cap- 
tured. Horizontal transport would be more important in a 3D field, leading to somewhat higher biases (Chambers, 
1997). Finally, only one wavelength, with no absorption, has been considered. Nevertheless, the results provide 
some insight into the effect of broken clouds on albedo bias. 

Section 2 briefly describes the two dimensional cloud fields inferred from LANDSAT. Section 3 describes the 
numerical radiation model as well as the radiative quantities studied in this work. Section 4 compares the cloud 
albedo biases based on the IPA and 2D results. Section 5 presents the major conclusions from this work. 

2. LANDSAT Inferred Horizontal Inhomogeneity 

Forty-five LANDSAT scenes of cloud fields over ocean were available for this study. The criteria that governed 
the selection of these particular scenes, and their characteristics, are described in detail in Harshvardhan et al. 
(1994) and Chambers et al. (1996). LANDSAT provides both better resolution and far more data volume than 
ground-based measurements of cloud properties. While the LANDSAT data may suffer from radiative smoothing 
effects at small scales (e.g., Marshak et al., 1995), test cases have shown that this does not affect results at the 
averaging scales used in this paper (Chambers et al., 1996; hereafter CWE). Note that seven of the scenes from 
Harshvardhan et al. (1994) were rejected because they had more than 10 percent of the pixels saturated in radiance. 
The final set of 45 scenes contains a majority that are subtropical oceanic boundary layer cloud, and therefore should 
not be interpreted in any way as a globally representative set. 

Based on the LANDSAT reflectance data, interpreted using IPA, the cloud optical depth, r, at 0.83 //m has been 
derived in each pixel as described in Harshvardhan et al. (1994). Using an empirical relation from Minnis et al. 
(1992), an estimate of the physical thickness, A z, of the cloud in each pixel can be made. The relation for A z was 
derived from simultaneous satellite and surface observations at a 30 km scale, but there is some evidence (CWE) 
that it is reasonable at the LANDSAT pixel scale as well. The extinction coefficient, /3, can then be calculated for 
each pixel and distributed through the cloud depth. For this study, the cloud is assumed to have uniform extinction 
with height in each pixel, i.e., (3 = r/Az. This extinction field is then discretized onto the computational grid, 
while maintaining total column optical depth at the LANDSAT value (see CWE). The cloud top height is held 
constant, while the cloud base varies according to the calculated Az. For inversion-capped boundary layer clouds 
over ocean, this is considered more physically realistic than a constant base height assumption. 

Conservative scattering (wo = 1) has been assumed throughout this study. A Mie phase function for a gamma 
distribution of water droplets with an effective radius, r e , of 10 /im and an effective variance of 0.1 is used and is 
assumed constant throughout the cloud. The mean cosine of the scattering phase function is 0.856 which is close 
to the value used by Cahalan et al. (1994b) for the Henyey-Greenstein phase function. The effect of atmospheric 
gases on the radiation within the cloud is neglected. While this may affect the initial retrieval of optical depth 
(Oreopoulos and Davies, 1996), it has little effect on horizontal radiative transfer characteristics and thus the 
conclusions of this paper (Chambers, 1997). 

From each LANDSAT scene, a series of 20 randomly placed horizontal scan line samples is taken. The orientation 
of clouds in the scenes is random enough that horizontal sampling is sufficient. Each scan line sample is 200 pixels 
(5.7 km) long. Sensitivity studies confirm this is a sufficient sample length for thin boundary layer clouds. The 
sampling is required because computing 2D transport in the entire scene is currently beyond the bounds of computer 
resources in terms of both time and memory. However PP-IPA bias can be calculated easily for the full scene and 
compared to the bias found from the sampled scans to test convergence. The correlation between the two is 
A a $ ampie = 0.98* A otfuii + 0.004, with a regression coefficient of 0.95 and a standard error of estimate of 0.01, 
indicating excellent convergence. What differences exist are traceable to a difference in mean properties between the 
sample and the full scene, and are accounted for since all results herein are examined in terms of cloud properties. 
A cyclical horizontal boundary condition is imposed in the 2D solution, with three logarithmic interpolating points 
added on the right side of most scans to blend the cloud field from edge to edge. In a few cases, additional 
interpolating points are added to keep the edge gradient no greater than the largest interior gradient of the scene. 

For the scan line samples taken from the LANDSAT scenes, mean cloud optical depth ranges from 0.3 to 
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73.3. Cloud fraction, A c ranges from 0.0 to 1.0. Aspect ratio (defined in Section 3c.) ranges from 0.0 to 2.73. 
Standard deviation of optical depth is between 0.08 and 45. The gamma function parameter (defined in Section 
3d.), measuring the width of the optical depth distribution, ranges from 0.04 to 310. 


3. Radiative Model 


a. SHDOM 

The Spherical Harmonic Discrete Ordinate Method (SHDOM) of Evans is used here in both an independent 
pixel and two-dimensional solution mode. More details on this method are given in CWE. The algorithm is an 
earlier version of a 3D method described in Evans (1996). In brief, it uses both spherical harmonics and discrete 
ordinates to represent the radiance field during different parts of the solution algorithm. The spherical harmonics 
are employed for efficiently computing the source function including the scattering integral. The discrete ordinates 
are used to integrate the radiative transfer equation through the spatial grid. The solution method is to simply 
iterate between the source function and radiance field, akin to a successive order of scattering approach. 

The two-dimensional (# and z) cloud fields generated from the LANDSAT scenes are provided as input to the 
SHDOM code. A solution is obtained with L = 11 zenith terms and M = 11 azimuthal terms for the spherical 
harmonic truncation; and with iV^ — 12 zenith angles and N <f> = 24 azimuth angles in the discrete ordinate 
discretization. Using the 6-M method (Wiscombe, 1977) for the phase function, this truncation has been found to 
provide good accuracy (mean error about 0.5 percent) for fluxes, upon which the results of this paper depend. A 
Monte Carlo method has been used to verify the results for two typical broken cloud scan line samples, showing 
agreement within the uncertainty of the Monte Carlo solution using 20 million photons. Along with prior work by 
one author (Evans), this is considered sufficient verification of the SHDOM numerical model. 

For each scan line sample of each scene, solutions are obtained for solar zenith angles, 0o> of 0°, 49° and 63° 
(^o = cos 0 O =1, 0.65, 0.45); and with both the IPA and 2D solution options. 


b. Albedo and Albedo Bias 


The spectral albedo, a(x,6a), is given by 


a(x, 9 q ) 


f T (M q) 

PqFo 


(1) 


where F© is the solar flux constant at the appropriate wavelength and F t is the upward spectral flux. All results 
in this paper are restricted to monochromatic calculations at 0.83 /im. F T is computed by SHDOM as a function 
of x for the 2D scene using both a 2D method and the IPA. To obtain the scan line sample average IPA and 2D 
albedos, then, a(x) is averaged over cloudy pixels in the 200 pixel sample. The PP albedo is retrieved from the 
average cloud optical depth of the scan (i.e., set r c (x) = f c for all x corresponding to cloudy pixels), using the 
results for albedo versus optical depth as shown in Fig. 1. Averages at larger scales are obtained analogously. 

The IPA-2D albedo bias for cloudy pixels only is defined as 


Aa/PA = (a ipA — ®2d) = 


^p%( a IPA ~ <* 2 d) 


while the PP-IPA bias for cloud areas is 


Aapp = (app — oc i pa) = ot PP (r c ) — 


Ep %ajPA 

Ep % 


(2) 

( 3 ) 


where = 0 for r — 0, and $ p = 1 for r > 0. The PP-2D bias is therefore the sum of (2) and (3). Albedo bias for 
the whole field, applicable to some General Circulation Models (GCM), is found by multiplying either definition by 
the cloud fraction. To obtain the relative bias, (2) is divided by aiPA\ (3) by app. These definitions are consistent 
with the work of Cahalan et al. (1994a,b). As in that work, relative bias is given here in percent; while absolute 
bias is not. 
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c. Aspect Ratio 

In order to classify the cloud fields in this study, a scan average cloud aspect ratio is computed according to 

CWE as: 

V_\ _ Az max 

h) 2Ax(p(Az) = 1/e) 

where p(Az) is the autocorrelation coefficient of the cloud thickness (i.e., the degree to which A z(x) and A z(x + dx) 
are correlated). This definition measures the ratio of the largest vertical cloud extent (A z max ) to a horizontal width 
which approximates the largest feature width at half-maximum for the sample (see CWE for more details and an 
example). Scan line samples which are completely clear over the 5.7 km sample are assigned an aspect ratio of 
zero. The aspect ratio for completely overcast scan line samples is computed from the largest vertical cloud extent 
divided by the full 5.7 km sample width. Other definitions are certainly possible, but this one has been found to 
be useful. 

d. Gamma Function Parameter, v 

The gamma function has been proposed by Barker et al. (1996) as an analytic approximation for the distribution 
of optical depths in an actual cloud field. Such an analytic form can enable an approximate IPA solution to be 
computed for little more than the cost of a PP result. The shape of the distribution is controlled by the parameter 
v , which is given by 



where r c is the mean and a T the standard deviation of the distribution of optical depths. Small values of v produce 
optical depths distributed in a monotonically decreasing fashion, characteristic of broken cloud fields. Values of v 
larger than one produce distributions with a more Gaussian shape and a well-defined mean, characteristic of overcast 
clouds. This parameter can also be viewed as a measure of the relative variability of a scene (u < 1 => cr T > f c => 
highly variable). 

4. Bias of PP and IPA Cloud Albedos 

At the pixel level the variability between IPA and 2D albedo is large, due to horizontal transport displacing the 
peaks and valleys in the radiative flux, and thus masking albedo biases. The results are therefore first averaged 
over the cloudy region of a scan line sample. The albedo bias over such small horizontal scales is still large and 
uncorrelated to cloud field properties (as found by Cahalan et al. (1994b)), due again to horizontal transport. The 
behavior of the bias becomes reasonable only on the scale of a LANDSAT scene (58 by 58 km) or even larger. 
The mean cloud albedos for the 20 scan line samples in each LANDSAT scene are therefore averaged together in 
several ways to provide a mesoscale or GCM scale view. The results are first presented on a scene average basis, 
in which results for the 20 scan line samples in each LANDSAT scene are averaged. The second averaging method 
classifies each scan line sample in terms of a cloud parameter, then averages all the samples from any LANDSAT 
scene in that cloud parameter bin. Results using the first averaging method still show considerable scatter when 
plotted against any cloud parameter, due to within-scene variability of the cloud field and sensitivity to variations 
in the other parameters within the cloud scene. The second approach allows the variability to be reduced using 
conditional sampling as a function of cloud parameters (cloud fraction A c , optical depth t c , aspect ratio, standard 
deviation a T , or gamma parameter u) and provides a clearer picture of the trend of albedo bias as a function of 
these cloud physical properties. The second approach has the same statistical distribution as the scene average 
results, but provides better sampling for uncertainty estimation. Finally, averaging is done in three broad cloud 
classes to generate mesoscale and climatological scale statistics. 

a. Scene-Average Results 

Albedo biases averaged for the 20 scan line samples in each scene represent the bias at a (58 km) 2 scale, applicable 
to GCM modeling. Results are presented as a function of the standard deviation of optical depth, a T , of the scene. 
This is the single parameter which best correlates with the PP-IPA bias (aspect ratio best correlates with the 
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IPA-2D bias). In terms of the other cloud physical parameters, the variation is similar to that shown in Sect. 4. b., 
generally falling within the uncertainty (1-sigma) of those results. As shown in Fig. 2a, the PP-IPA bias generally 
increases with a T , i.e., with the variability of the optical depth field, while the IPA-2D bias is an order of magnitude 
or more smaller. The PP-IPA bias is strictly positive, due to the shape of the curve in Fig. 1, while at overhead 
sun the IPA-2D bias is almost evenly distributed among positive and negative values: positive (IPA > 2D) for the 
broken cloud scenes in which leakage of photons out cloud sides and to the surface dominates; negative (IPA < 
2D) for more overcast ones in which “filling in” (i.e., reflected radiance from cloud sides that emerges at the top of 
a cloud-free area) between closely spaced clouds or cloud cells dominates. 

At the higher sun angles shown in Figs. 2b and 2c the results are similar, with somewhat more scatter - 
particularly in the IPA-2D bias - as a r increases. The PP-IPA bias remains strictly positive while the IPA-2D bias 
now tends to more negative values, particularly for the more variable scenes. This is due to the fact that at slant 
sun the cloud sides intercept and reflect some of the radiation that impinges on cloud-free regions at the top of the 
field, making the 2D albedo larger than the IPA albedo (filling-in effect). At 9q — 63°, this effect is overwhelmed 
for some of the more uniform scenes by shadowing of cloud sides in 2D, resulting in a larger IPA albedo and a 
positive bias. 

The importance of horizontal radiative transport to albedo bias can be estimated from the size of the IPA-2D 
bias relative to the PP-2D bias. For the scene averages, the IPA-2D bias accounts for 20% Root- Mean- Square 
(RMS) of the albedo bias, indicating horizontal transport is generally a small - but not necessarily negligible - 
factor in the albedo bias. 

b. Parameter Space Results 

Albedo biases averaged over the 5.7 km scan line samples are binned and averaged to examine the sensitivity of 
the bias to cloud properties. Results are shown in Figs. 3a to 3e in terms of five cloud parameters: cloud fraction, 
mean cloud optical depth, cloud aspect ratio, standard deviation of optical depth, and gamma function parameter, 
v. The cloud fraction bins are those planned for use in the CERES experiment ( Baum et a/., 1995), while bins 
for the other cloud properties were chosen after examining the data. The trends shown are not very sensitive to 
bin boundaries. Note that the PP-IPA bias can theoretically depend only on the frequency distribution of optical 
depth (i.e., r, <r r and z/). The apparent dependence on cloud fraction and aspect ratio comes about because A c 
and aspect ratio are correlated with f, a T or v in this set of cloud scenes. (For example, small A c tends to mean 
thin clouds in this dataset.) 

Two sets of error bars are given in these figures, for the 63° solar zenith angle results where the error is generally 
largest. The error bar through the data points is an estimate of the error in the mean value. It is generated 
by splitting the cloud scenes into two completely uncorrelated sets (scenes from different years) and taking the 
difference in the mean value found for each set. This error is quite small, except in bins where there are very few 
(about 20 or fewer) samples, suggesting good confidence in the mean bias value. The second error bar, slightly 
offset from the data points, represents the cloud fraction weighted 1-sigma variability of the scan line samples in 
that bin. 

The variation with cloud fraction is shown in Fig. 3a. The PP-IPA bias is near zero for small cloud fractions, 
which in these LANDSAT scenes correspond to thin clouds at the linear limit of the retrieval curve in Fig. 1. The 
bias increases for intermediate cloud fractions, then falls again in the overcast limit. This behavior is due to the 
correlation with cloud field variability in this dataset, v is generally larger, implying less variability, when the cloud 
fraction is 1.0, compared to the broken cloud fields in this set. The IPA-2D bias has a similar trend, but with a 
much smaller magnitude and the opposite sign. The trends are similar at all three solar zenith angles, though the 
bias tends to be smaller at overhead sun (0o — 0°). At slant sun, the IPA-2D bias again is consistent with added 
energy incident into the sides of broken clouds leading to a larger 2D albedo. 

In terms of optical depth, Fig. 3b, the PP-IPA bias starts out very small as it should in the linear limit. If the 
optical depth ranged high enough, this bias should drop again at large r c , as the retrieval curve (Fig. 1) becomes 
nearly linear again around r = 50 (unless the variability also increased with r c ). Since in this sample the maximum 
optical depth is below 30, the PP-IPA bias remains quite large. The IPA-2D bias, on the other hand, peaks below 
r = 5 before dropping again; most obviously at 0q = 63°. 



6 


The results as a function of aspect ratio are shown in Fig. 3c. The PP-IPA bias increases almost linearly with 
aspect ratio for all solar zenith angles. The IPA-2D bias has an even more linear trend. Aspect ratio is in fact 
the parameter that best correlates the IPA-2D scan line sample bias (r = 0.56 to 0.79 as 6$ increases), suggesting 
that cloud geometry is somewhat more important than the features of the optical depth distribution in controlling 
horizontal radiative transport. 

In terms of the standard deviation of optical depth, <j t , results are shown in Fig. 3d. The PP-IPA bias in this 
case also increases quite linearly with a T , at all sun angles. This is the parameter that best correlates the PP-IPA 
scan line sample bias, with a regression coefficient between 0.71 and 0.88 depending on sun angle. The IPA-2D 
bias shows a similar behavior but with an increasingly negative bias as the sun drops from overhead for the more 
variable scenes due to energy intercepted by cloud sides. 

Finally, the trends in terms of gamma function parameter, v , can be observed from Fig. 3e. The biases are nearly 
constant for large v . As v decreases below one, implying a highly variable cloud field, the biases increase rapidly. 
In fact, based on the magnitude of the largest biases in this figure, v does the best job of grouping scans with large 
Aa. This suggests that v is a powerful parameter in determining the influence of horizontal cloud inhomogeneity, 
and should be considered in cloud field classification studies. 

Overall, as found by Cahalan et al. (1994a), the albedo bias is less sensitive to mean properties (A c , r) than it 
is to properties of the distribution of r (<r r , z/, and aspect ratio). 

All the parameter space results show the IPA-2D bias to be smaller - sometimes by orders of magnitude - than 
the PP-IPA bias, except at small optical depth or cloud fraction where the cloud field approaches the linear limit; 
and at large aspect ratio or small v where the IPA-2D bias can be half the PP-IPA bias at slant sun. 

c. Cloud Amount Class Results 

To examine larger scale results, the LANDSAT scenes are now composited together in several ways. Results for 
the whole dataset are summarized in Tables la and lb, in terms of three cloud amount classes roughly representing 
scattered, broken, and overcast cloud fields. 

Table la summarizes the mean and 1-sigma variability of the scene average ((58 km) 2 ) properties for all scenes 
in a cloud class. (Similar results for spherical albedo bias in Barker et al. (1996) use the same initial dataset 
but slightly different categorization and scene rejection criteria.) Such results are most applicable to a GCM 
or mesoscale model. The albedo bias and the variability in that bias are found to be largest for the smallest 
cloud fractions. This is consistent with the amount of variability in each of the cloud classes, as measured by 
<7 T and za Cahalan et al, (1994a), in measurements during FIRE [First ISCCP (International Satellite Cloud 
Climatology Project) Regional Experiment] off the coast of California, found the largest variability in overcast 
clouds. In measurements taken during ASTEX (Atlantic Stratocumulus Experiment), on the other hand, Cahalan 
et al. (1995) found more variability in broken clouds. More data, obtained from a variety of instruments, may be 
required to better understand cloud variability on a global basis. 

Much of the work of Cahalan et al. (1994a,b) was done using data collected over a period of several days or 
weeks. To better compare with long-term average data of this sort, Table lb groups the LANDSAT scenes in each 
cloud amount class into composite “scattered”, “broken”, and “overcast” cloud scenes, and presents the albedo 
bias for each (akin to (225 km) 2 area averages). It also gives results for a single cloud scene compositing all 45 
LANDSAT scenes (“all scenes”, (400 km) 2 average). The results obtained by Cahalan et al. (1994a) correspond to 
mostly cloudy data between the “overcast” and “all scenes” averages, and at 63° solar zenith angle these results 
agree well with their 15% bias. The 8% bias found in Marshak et al. (1995) is for overcast clouds at an averaging 
scale between Table la and Table lb and also agrees well with the present results. 

Examination of the relative magnitude of the IPA-2D versus PP-2D bias shows that the effect of horizontal 
radiative transport on albedo bias is largest for scattered clouds and smallest for overcast clouds. IPA-2D bias 
accounts for 32% (RMS over the scenes) of the bias for the scattered cloud class, 11% RMS for broken cloud fields, 
and 7% RMS in the overcast clouds. The overcast result is identical to the relative bias contributions found in 
Cahalan et al. (1994b). 

Recall that this dataset is not meant to represent the global distribution of clouds; and these bias numbers should 
therefore not be viewed in any way as global averages. The data should be representative of marine boundary layer 
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clouds, however. 

5. Conclusions 

The albedo bias of the independent pixel approximation (IPA-2D) and plane parallel method (PP-IPA) for 
oceanic boundary layer clouds have been assessed using realistic cloud fields from 45 LANDSAT scenes including 
broken cloudiness, and a two-dimensional SHDOM radiative transport method. The modeling includes conservative 
scattering, ignores gas absorption, and uses a black surface. Tests indicate that these assumptions have only a minor 
effect on the results shown. 

PP-IPA cloud albedo biases up to near 0.30 were found in this study. They are strictly positive due to the shape 
of the albedo versus optical depth relationship and are best correlated by the standard deviation of the optical 
depth distribution, <r T . IPA-2D bias, on the other hand, ranges from very small positive values to near -0.12 for 
extremely variable cloud fields. Cloud aspect ratio, as defined in this study, is found to best correlate this bias. 
While the IPA-2D bias is generally much smaller than the PP-IPA bias, it is not always negligible. Overall, it 
accounts for 20% RMS of the bias. Negative IPA-2D bias occurs due to reflection from the cloud sides in the 
2D solution. The dependence of albedo bias on five cloud parameters (cloud fraction, mean cloud optical depth, 
cloud aspect ratio, standard deviation of optical depth, and gamma function shape parameter) is also examined. 
The most useful parameters for predicting bias due to cloud field horizontal inhomogeneities are found to be the 
standard deviation of the optical depth distribution, the gamma function parameter and the cloud aspect ratio. 
Finally, the bias is examined at larger mesoscale and climatological scales, and found to be in good agreement with 
earlier work on overcast and mostly cloudy cases. New results are given here for broken and scattered cloud fields 
which show the albedo bias can be even larger in these conditions, depending on the within cloud variability of the 
cloud field. 
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